
library(pacman)
pacman::p_load("tidyverse","lavaan","kableExtra","xtable","ggcorrplot","gridExtra","riclpmr")

load("Data/gesis_cleaned.RData")

#calculate
nrow(na.omit(dplyr::select(gesis,open15_latent,open17_latent,open19_latent,lr19, lr17, lr15)))
nrow(na.omit(dplyr::select(gesis,con15_latent,con17_latent,con19_latent,lr19, lr17, lr15)))

nrow(na.omit(dplyr::select(gesis,open15_latent,open19_latent,antiEU14, antiEU18)))
nrow(na.omit(dplyr::select(gesis,con15_latent,con19_latent,antiEU14, antiEU18)))

#rename dataset
all <- gesis_long

#Figure 3 -------------
#set model
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1 + variable
  x3 + y3 ~ x2 + y2 + variable
  x3 ~~ x3
  x2 ~~ x2
  x1 ~~ x1
  y3 ~~ y3
  y2 ~~ y2
  y1 ~~ y1
  y3 ~~ x3
  y2 ~~ x2
  y1 ~~ x1
  y1 ~ variable
'


CLPM_c <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ a*x1 + b*y1 + p*variable
  x3 + y3 ~ a*x2 + b*y2 + p*variable
  x3 ~~ x3
  x2 ~~ x2
  x1 ~~ x1
  y3 ~~ y3
  y2 ~~ y2
  y1 ~~ y1
  y3 ~~ x3
  y2 ~~ x2
  y1 ~~ x1
  y1 ~ variable

'

#select data and prepare variables for left-right self-placement model
xvars <- colnames(all)[grep("lr", x = colnames(all))][c(1,2,3)]
yvars <- colnames(all)[grep("p", x = colnames(all))][c(8,9,10)]
all %>% select(xvars, yvars, z000001a, variable) -> gss_short1

names(gss_short1)[1:3] <- paste0("x", 1:3)
names(gss_short1)[4:6] <- paste0("y", 1:3)
gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$x3<- scale(gss_short1$x3)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$y3<- scale(gss_short1$y3)

#run models
lr_clpm <- sem(CLPM_u, data = (gss_short1),missing="fiml",fixed.x=F, cluster="z000001a")
lr_clpm_con <- sem(CLPM_c, data = (gss_short1),missing="fiml",fixed.x=F, cluster="z000001a")

#compare unconstrained vs. contrained model
anova(lr_clpm,lr_clpm_con)

#set model for EU
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1 + variable
  x2 ~~ x2
  x1 ~~ x1
  y2 ~~ y2
  y1 ~~ y1
  y2 ~~ x2
  y1 ~~ x1
  y1 ~ variable
'

#select data and prepare variables for EU model
xvars <- colnames(all)[grep("antiEU", x = colnames(all))][c(1,2)]
yvars <- colnames(all)[grep("p", x = colnames(all))][c(9,10)]
all %>% select(xvars, yvars, z000001a, variable) -> gss_short1
names(gss_short1)[1:2] <- paste0("x", 1:2)
names(gss_short1)[3:4] <- paste0("y", 1:2)
gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
eu_clpm <- sem(CLPM_u, data = (gss_short1),missing="fiml",fixed.x=F,cluster = "z000001a")

#create figure
auth_results <- rbind(parameterestimates(lr_clpm,level = .90),parameterestimates(eu_clpm,level = .90)) %>% filter(lhs %in% c("x2","x3","y2","y3") & rhs %in% c("x2","x1","y2","y1") & op == "~") %>% slice(c(2,3,6,7,10,11))
auth_results$lhs <- c("DV: Policy Preference","DV: Personality")
auth_results$rhs <- c("IV: Personality","IV: Policy Preference")
auth_results$label <- c(rep("Ideology",4),rep("EU",2))
auth_results$lhs <- factor(auth_results$lhs,levels=c("DV: Policy Preference","DV: Personality"))
auth_results$wave <- c(rep("Wave 2 DV~Wave 1 IV",2),rep("Wave 3 DV~Wave 2 IV",2),rep("Wave 2 DV~Wave 1 IV",2))
ggplot(auth_results,aes(x=label,y=est,shape=wave))+geom_pointrange(aes(ymin=ci.lower,ymax=ci.upper),position=position_dodge(width=.5))+facet_wrap(lhs~rhs)+geom_hline(yintercept = 0,lty=2)+theme_bw(base_size = 15)+labs(shape="Model")+ylab("Standardized Coefficient")+xlab("Policy Variable") + theme(legend.position="bottom")
ggsave("Results/xlaggedplotgesis.pdf",width=8,height=6)

#SI 2.6 - Results belonging to Figure 3 in the main text --------
# For l-r
parameterEstimates(lr_clpm) %>% filter((str_detect(lhs, "^x") & str_detect(rhs, "^y")) | (str_detect(lhs, "^y") & str_detect(rhs, "^x")) ) %>% select(1,3,4,5,7) %>% slice(c(1:4)) -> euth_param
euth_param$pvalue <- euth_param$pvalue/2
fitindicies <-
  rbind(c(
    "RMSEA",
    round(fitmeasures(lr_clpm)[names(fitmeasures(lr_clpm)) == "rmsea"], 2),
    paste0("CI:",round(fitmeasures(lr_clpm)[names(fitmeasures(lr_clpm)) ==
                                              "rmsea.ci.lower"], 2),"-",
           round(fitmeasures(lr_clpm)[names(fitmeasures(lr_clpm)) ==
                                        "rmsea.ci.upper"], 2)),NA, NA
  ),
  c("CFI", round(fitmeasures(lr_clpm)[names(fitmeasures(lr_clpm)) ==
                                        "cfi"], 2), rep(NA, 3)),
  c("SRMR", round(fitmeasures(lr_clpm)[names(fitmeasures(lr_clpm)) ==
                                         "srmr"], 2), rep(NA, 3)))
colnames(fitindicies) <- colnames(euth_param)
euth_tab <- rbind(euth_param[1:4,],fitindicies)
euth_tab[1:4,3:5] <- mapply(function(x)as.character(formatC(round(as.numeric(x),2),2,format="f")), euth_tab[1:4,3:5])

euth_tab$lhs <- c("Politics", "Personality","Politics", "Personality","RMSEA","CFI","SRMR")
euth_tab$rhs <- c("$Personality_{t-1}$", "$Politics_{t-1}$",  "$Personality_{t-1}$","$Politics_{t-1}$", euth_tab$rhs[5:7])
colnames(euth_tab) <- c("Outcome Measure","Predictor","$\\beta$","SE","p-value")
kbl(euth_tab,booktabs = T, digits=2, caption="Cross-lagged coefficients between Left-Right Self-Placement and Closed Personality", format = "latex", escape = F) %>%
  kable_classic(full_width = F) %>%
  kableExtra::group_rows(start_row = 5,end_row = 7,group_label = "Fit Indices") %>%
  kableExtra::group_rows(start_row = 1,end_row = 2,group_label = "Wave 2") %>%
  kableExtra::group_rows(start_row = 3,end_row = 4,group_label = "Wave 3") %>% kable_styling(latex_options = c( "hold_position")) %>%  kable_styling(latex_options = c( "hold_position")) ->  lr_tab
lr_tab %>% save_kable("Results/gesis_lr.tex")

# For EU
parameterEstimates(eu_clpm) %>% filter((str_detect(lhs, "^x") & str_detect(rhs, "^y")) | (str_detect(lhs, "^y") & str_detect(rhs, "^x")) ) %>% select(c(1,3:5,7)) -> libcon_param
eu_clpm <- data.frame(mapply(function(x)as.character(formatC(round(x,2),2,format="f")), libcon_param[1:2,3:5])) -> eu_tab

eu_clpm$lhs <- c("Politics", "Personality")
eu_clpm$rhs <- c("$Personality_{t-1}$", "$Politics_{t-1}$")
eu_clpm %>% relocate(4:5, .before=est) %>% mutate(pvalue=as.numeric(as.character(pvalue))/2) -> eu_clpm
colnames(eu_clpm) <- c("Outcome Measure","Predictor","$\\beta$","SE","p-value")
kbl(eu_clpm,booktabs = T, digits=2, caption="Cross-lagged coefficients between EU Attitudes and Closed Personality", format = "latex", escape = F) %>%
  kable_classic(full_width = F) %>% kable_styling(latex_options = c( "hold_position")) -> eu_tab
eu_tab %>% save_kable("Results/gesis_eu.tex")


#SI 2.3 - Psychometric properties------------

wave1psych <-
  list(
    psych::alpha(data.frame(gesis$open15_1, gesis$open15_2))[[1]],
    psych::alpha(data.frame(gesis$con15_1, gesis$con15_2))[[1]],
    psych::alpha(data.frame(gesis$ext15_1, gesis$ext15_2))[[1]],
    psych::alpha(data.frame(gesis$agre15_1, gesis$agre15_2))[[1]],
    psych::alpha(data.frame(gesis$neu15_1, gesis$neu15_2))[[1]]
  )

wave2psych <-
  list(
    psych::alpha(data.frame(gesis$open17_1, gesis$open17_2))[[1]],
    psych::alpha(data.frame(gesis$con17_1, gesis$con17_2))[[1]],
    psych::alpha(data.frame(gesis$ext17_1, gesis$ext17_2))[[1]],
    psych::alpha(data.frame(gesis$agre17_1, gesis$agre17_2))[[1]],
    psych::alpha(data.frame(gesis$neu17_1, gesis$neu17_2))[[1]]
  )

wave3psych <-
  list(
    psych::alpha(data.frame(gesis$open19_1, gesis$open19_2))[[1]],
    psych::alpha(data.frame(gesis$con19_1, gesis$con19_2))[[1]],
    psych::alpha(data.frame(gesis$ext19_1, gesis$ext19_2))[[1]],
    psych::alpha(data.frame(gesis$agre19_1, gesis$agre19_2))[[1]],
    psych::alpha(data.frame(gesis$neu19_1, gesis$neu19_2))[[1]]
  )

wave1 <- do.call("rbind",wave1psych)[,c(2,4)]
wave2 <- do.call("rbind",wave2psych)[,c(2,4)]
wave3 <- do.call("rbind",wave3psych)[,c(2,4)]

alphas <- rbind(data.frame(wave=2015, trait= NA ,wave1),
                data.frame(wave=2017, trait= NA ,wave2),
                data.frame(wave=2019, trait= NA ,wave3))
alphas$trait <- rep(c("Openness","Conscientiousness","Extraversion","Agreeableness","Neuroticism"),3)
colnames(alphas) <- c("Wave","Trait","Alpha","Average r")
print(xtable(alphas, caption="Reliability of Personality Measures by Wave", digits=c(0,0,0, 2,2), align=c("llccc")),include.rownames=F,caption.placement="top",file="Results/appendix_gesis_alphas.tex")

#SI 2.4 - descriptives: correlation between variables----------------
w1 <- mapply(zero1,with(gesis,data.frame(lr15, antiEU14, open15_latent, con15_latent, ext15_latent, agre15_latent, neu15_latent)))
w2 <- mapply(zero1,with(gesis,data.frame(lr17, open17_latent, con17_latent, ext17_latent, agre17_latent, neu17_latent)))
w3 <- mapply(zero1,with(gesis,data.frame(lr19, antiEU18, open19_latent, con19_latent, ext19_latent, agre19_latent, neu19_latent)))

reshape2::melt(data.frame(w1)) -> w1a #ignore warning
w1a %>% group_by(variable) %>% summarise(mean=mean(value,na.rm=T),sd=sd(value,na.rm=T),N=sum(value,na.rm=T)) -> w1a
reshape2::melt(data.frame(w2)) -> w2a #ignore warning
w2a %>% group_by(variable) %>% summarise(mean=mean(value,na.rm=T),sd=sd(value,na.rm=T),N=sum(value,na.rm=T)) -> w2a
reshape2::melt(data.frame(w3)) -> w3a #ignore warning
w3a %>% group_by(variable) %>% summarise(mean=mean(value,na.rm=T),sd=sd(value,na.rm=T),N=sum(value,na.rm=T)) -> w3a

alldescr <- rbind(w1a,w2a,w3a)
alldescr <- data.frame(wave=c(2015, 2015, 2015, 2015, 2015, 2015, 2015, 2017, 2017, 2017, 2017, 2017, 2017, 2019, 2019, 2019, 2019, 2019, 2019, 2019),variable=c("Liberal-Conservative", "EU", "Openness","Conscientiousness","Extraversion","Agreeableness","Neuroticism", "Liberal-Conservative", "Openness","Conscientiousness","Extraversion","Agreeableness","Neuroticism", "Liberal-Conservative", "EU", "Openness","Conscientiousness","Extraversion","Agreeableness","Neuroticism"),alldescr[,-1])
print(xtable(alldescr, caption="Descriptive Statistics of Measures by Wave", digits=c(0,0,0, 2,2,0),align=c("llcccc")), label = "tab:gesis_descriptives",  include.rownames=F,caption.placement="top", file="Results/appendix_gesis_descri.tex")


#SI 2.5 - #descriptives: correlation between variables----------------
w1cor <- cor(w1,use="pairwise.complete.obs")
rownames(w1cor) <- c("Left-right", "EU", "Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism")
colnames(w1cor) <- c("Left-right", "EU", "Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism")
a<-ggcorrplot(w1cor, type="lower", show.diag = TRUE, lab=TRUE, lab_size = 2 , title="A. Measures from 2015")

w2cor <- cor(w2,use="pairwise.complete.obs")
rownames(w2cor) <- c("Left-right", "Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism")
colnames(w2cor) <- c("Left-right", "Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism")
b<-ggcorrplot(w2cor, type="lower", show.diag = TRUE, lab=TRUE, lab_size = 2, title="B. Measures from 2017")

w3cor <- cor(w3,use="pairwise.complete.obs")
rownames(w3cor) <- c("Left-right", "EU", "Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism")
colnames(w3cor) <- c("Left-right", "EU", "Openness", "Conscientiousness", "Extraversion", "Agreeableness", "Neuroticism")
c<-ggcorrplot(w3cor, type="lower", show.diag = TRUE, lab=TRUE, lab_size = 2, title="C. Measures from 2019")

corplot <- grid.arrange(a,b, c)
ggsave(corplot,file="Results/gesis_correlations.pdf",width=20,height=20)

#SI 2.7 - Crossed-lagged coefficients for RI-CLPM models-------------

#select and prepare data
riclpmtext <- riclpm_text(var_groups = list(x=c("x1","x2","x3"),y=c("y1","y2","y3")))
xvars <- colnames(all)[grep("lr", x = colnames(all))][c(1,2,3)]
yvars <- colnames(all)[grep("p", x = colnames(all))][c(8,9,10)]
all %>% select(xvars, yvars,id,female,starts_with("education"),starts_with("income"),age,variable) -> gss_short1
names(gss_short1)[1:3] <- paste0("x", 1:3)
names(gss_short1)[4:6] <- paste0("y", 1:3)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$x3<- scale(gss_short1$x3)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$y3<- scale(gss_short1$y3)
gss_short1$id <- 1:nrow(gss_short1)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$x3<- scale(gss_short1$x3)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$y3<- scale(gss_short1$y3)
gss_short1$id <- 1:nrow(gss_short1)

#run RI-CLPM model
lr_tab <- riclpmr::lavriclpm(riclpmtext,data = gss_short1,missing="ML")

#extract coefficients and create table
parameterEstimates(lr_tab) %>% filter(label %in% c("x_y","y_x"))%>%  slice(1,3) -> riclpmliss
riclpmliss$pvalue <- riclpmliss$pvalue/2
riclpmliss$dv <- c("Politics","Personality")
riclpmliss$iv <- c("$Personality_{t-1}$", "$Politics_{t-1}$")
riclpmliss %>% select(dv,iv,est,se,pvalue) -> riclpmliss
colnames(riclpmliss) <- c("Outcome Measure","Predictor","$\\beta$","SE","p-value")
kbl(riclpmliss,booktabs = T, digits=2, caption="Cross-lagged coefficients from RI-CLPM models in GESIS model", format = "latex", escape = F) %>%
  kable_classic(full_width = F) %>%  kable_styling(latex_options = c( "hold_position")) %>%
  save_kable("Results/gesisRICLPM.tex")

#SI 2.8 - Results for Openness separately--------------------------
#set model
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1
  x3 + y3 ~ x2 + y2
  x3 ~~ x3
  x2 ~~ x2
  x1 ~~ x1
  y3 ~~ y3
  y2 ~~ y2
  y1 ~~ y1
  y3 ~~ x3
  y2 ~~ x2
  y1 ~~ x1
'

#select data and run model for Left-right placement
xvars <- colnames(gesis)[grep("lr", x = colnames(gesis))][c(1,2,3)]
yvars <- colnames(gesis)[grep("open", x = colnames(gesis))][c(3,6,9)]
gesis %>% select(xvars, yvars) -> gss_short1
names(gss_short1)[1:3] <- paste0("x", 1:3)
names(gss_short1)[4:6] <- paste0("y", 1:3)
gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$x3<- scale(gss_short1$x3)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$y3<- scale(gss_short1$y3)
lr_clpm <- sem(CLPM_u, data = (gss_short1),missing="ML",fixed.x=F)

#set model for EU
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1
  x2 ~~ x2
  x1 ~~ x1
  y2 ~~ y2
  y1 ~~ y1
  y2 ~~ x2
  y1 ~~ x1
'
#select data and run model for EU model
xvars <- colnames(gesis)[grep("antiEU", x = colnames(gesis))][c(1,2)]
yvars <- colnames(gesis)[grep("open", x = colnames(gesis))][c(6,9)]
gesis %>% select(xvars, yvars) -> gss_short1
names(gss_short1)[1:2] <- paste0("x", 1:2)
names(gss_short1)[3:4] <- paste0("y", 1:2)
gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
eu_clpm <- sem(CLPM_u, data = (gss_short1),missing="fiml",fixed.x=F)

#create plot
auth_results <- rbind(parameterestimates(lr_clpm,level = .90),parameterestimates(eu_clpm,level = .90)) %>% filter(lhs %in% c("x2","x3","y2","y3") & rhs %in% c("x2","x1","y2","y1") & op == "~") %>% slice(c(2,3,6,7,10,11))
auth_results$lhs <- c("DV: Policy Preference","DV: Personality")
auth_results$rhs <- c("IV: Personality","IV: Policy Preference")
auth_results$label <- c(rep("Ideology",4),rep("EU",2))
auth_results$lhs <- factor(auth_results$lhs,levels=c("DV: Policy Preference","DV: Personality"))
auth_results$wave <- c(rep("Wave 2 DV~Wave 1 IV",2),rep("Wave 3 DV~Wave 2 IV",2),rep("Wave 2 DV~Wave 1 IV",2))
ggplot(auth_results,aes(x=label,y=est,shape=wave))+geom_pointrange(aes(ymin=ci.lower,ymax=ci.upper),position=position_dodge(width=.5))+facet_wrap(lhs~rhs)+geom_hline(yintercept = 0,lty=2)+theme_bw(base_size = 15)+labs(shape="Model")+ylab("Standardized Coefficient")+xlab("Policy Variable") + theme(legend.position="bottom")
ggsave("Results/xlaggedplotgesis_openness.pdf",width=8,height=6)


#SI 2.8 - Results for Conscientiousness separately--------------------------
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1
  x3 + y3 ~ x2 + y2
  x3 ~~ x3
  x2 ~~ x2
  x1 ~~ x1
  y3 ~~ y3
  y2 ~~ y2
  y1 ~~ y1
  y3 ~~ x3
  y2 ~~ x2
  y1 ~~ x1
'

#Left-right
xvars <- colnames(gesis)[grep("lr", x = colnames(gesis))][c(1,2,3)]
yvars <- colnames(gesis)[grep("con", x = colnames(gesis))][c(3,6,9)]
gesis %>% select(xvars, yvars) -> gss_short1
names(gss_short1)[1:3] <- paste0("x", 1:3)
names(gss_short1)[4:6] <- paste0("y", 1:3)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$x3<- scale(gss_short1$x3)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$y3<- scale(gss_short1$y3)

lr_clpm <- sem(CLPM_u, data = (gss_short1),missing="ML",fixed.x=F)

#EU
xvars <- colnames(gesis)[grep("antiEU", x = colnames(gesis))][c(1,2)]
yvars <- colnames(gesis)[grep("con", x = colnames(gesis))][c(6,9)]
gesis %>% select(xvars, yvars) -> gss_short1
names(gss_short1)[1:2] <- paste0("x", 1:2)
names(gss_short1)[3:4] <- paste0("y", 1:2)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
#gss_short1$id <- 1:nrow(gss_short1)

CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1
  x2 ~~ x2
  x1 ~~ x1
  y2 ~~ y2
  y1 ~~ y1
  y2 ~~ x2
  y1 ~~ x1
'
eu_clpm <- sem(CLPM_u, data = (gss_short1),missing="fiml",fixed.x=F)

auth_results <- rbind(parameterestimates(lr_clpm,level = .90),parameterestimates(eu_clpm,level = .90)) %>% filter(lhs %in% c("x2","x3","y2","y3") & rhs %in% c("x2","x1","y2","y1") & op == "~") %>% slice(c(2,3,6,7,10,11))

auth_results$lhs <- c("DV: Policy Preference","DV: Personality")
auth_results$rhs <- c("IV: Personality","IV: Policy Preference")
auth_results$label <- c(rep("Ideology",4),rep("EU",2))
auth_results$lhs <- factor(auth_results$lhs,levels=c("DV: Policy Preference","DV: Personality"))
auth_results$wave <- c(rep("Wave 2 DV~Wave 1 IV",2),rep("Wave 3 DV~Wave 2 IV",2),rep("Wave 2 DV~Wave 1 IV",2))

ggplot(auth_results,aes(x=label,y=est,shape=wave))+geom_pointrange(aes(ymin=ci.lower,ymax=ci.upper),position=position_dodge(width=.5))+facet_wrap(lhs~rhs)+geom_hline(yintercept = 0,lty=2)+theme_bw(base_size = 15)+labs(shape="Model")+ylab("Standardized Coefficient")+xlab("Policy Variable") + theme(legend.position="bottom")
ggsave("Results/xlaggedplotgesis_conscientiousness.pdf",width=8,height=6)


#SI 2.9 - Exploratory results for Neuroticism, Neuroticism--------------------------
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1
  x3 + y3 ~ x2 + y2
  x3 ~~ x3
  x2 ~~ x2
  x1 ~~ x1
  y3 ~~ y3
  y2 ~~ y2
  y1 ~~ y1
  y3 ~~ x3
  y2 ~~ x2
  y1 ~~ x1
'

#Left-right
xvars <- colnames(gesis)[grep("lr", x = colnames(gesis))][c(1,2,3)]
yvars <- colnames(gesis)[grep("neu", x = colnames(gesis))][c(3,6,9)]
gesis %>% select(xvars, yvars) -> gss_short1
names(gss_short1)[1:3] <- paste0("x", 1:3)
names(gss_short1)[4:6] <- paste0("y", 1:3)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$x3<- scale(gss_short1$x3)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$y3<- scale(gss_short1$y3)

lr_clpm <- sem(CLPM_u, data = (gss_short1),missing="ML",fixed.x=F)

#EU
xvars <- colnames(gesis)[grep("antiEU", x = colnames(gesis))][c(1,2)]
yvars <- colnames(gesis)[grep("neu", x = colnames(gesis))][c(6,9)]
gesis %>% select(xvars, yvars) -> gss_short1
names(gss_short1)[1:2] <- paste0("x", 1:2)
names(gss_short1)[3:4] <- paste0("y", 1:2)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)

CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1
  x2 ~~ x2
  x1 ~~ x1
  y2 ~~ y2
  y1 ~~ y1
  y2 ~~ x2
  y1 ~~ x1
'
eu_clpm <- sem(CLPM_u, data = (gss_short1),missing="fiml",fixed.x=F)

auth_results <- rbind(parameterestimates(lr_clpm,level = .90),parameterestimates(eu_clpm,level = .90)) %>% filter(lhs %in% c("x2","x3","y2","y3") & rhs %in% c("x2","x1","y2","y1") & op == "~") %>% slice(c(2,3,6,7,10,11))

auth_results$lhs <- c("DV: Policy Preference","DV: Personality")
auth_results$rhs <- c("IV: Personality","IV: Policy Preference")
auth_results$label <- c(rep("Ideology",4),rep("EU",2))
auth_results$lhs <- factor(auth_results$lhs,levels=c("DV: Policy Preference","DV: Personality"))
auth_results$wave <- c(rep("Wave 2 DV~Wave 1 IV",2),rep("Wave 3 DV~Wave 2 IV",2),rep("Wave 2 DV~Wave 1 IV",2))

ggplot(auth_results,aes(x=label,y=est,shape=wave))+geom_pointrange(aes(ymin=ci.lower,ymax=ci.upper),position=position_dodge(width=.5))+facet_wrap(lhs~rhs)+geom_hline(yintercept = 0,lty=2)+theme_bw(base_size = 15)+labs(shape="Model")+ylab("Standardized Coefficient")+xlab("Policy Variable") + theme(legend.position="bottom")
ggsave("Results/xlaggedplotgesis_neuroticism.pdf",width=8,height=6)

#SI 2.9 - Exploratory results for Extraversion--------------------------
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1
  x3 + y3 ~ x2 + y2
  x3 ~~ x3
  x2 ~~ x2
  x1 ~~ x1
  y3 ~~ y3
  y2 ~~ y2
  y1 ~~ y1
  y3 ~~ x3
  y2 ~~ x2
  y1 ~~ x1
'

#Left-right
xvars <- colnames(gesis)[grep("lr", x = colnames(gesis))][c(1,2,3)]
yvars <- colnames(gesis)[grep("ext", x = colnames(gesis))][c(3,6,9)]
gesis %>% select(xvars, yvars) -> gss_short1
names(gss_short1)[1:3] <- paste0("x", 1:3)
names(gss_short1)[4:6] <- paste0("y", 1:3)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$x3<- scale(gss_short1$x3)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$y3<- scale(gss_short1$y3)

lr_clpm <- sem(CLPM_u, data = (gss_short1),missing="ML",fixed.x=F)

#EU
xvars <- colnames(gesis)[grep("antiEU", x = colnames(gesis))][c(1,2)]
yvars <- colnames(gesis)[grep("ext", x = colnames(gesis))][c(6,9)]
gesis %>% select(xvars, yvars) -> gss_short1
names(gss_short1)[1:2] <- paste0("x", 1:2)
names(gss_short1)[3:4] <- paste0("y", 1:2)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)

CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1
  x2 ~~ x2
  x1 ~~ x1
  y2 ~~ y2
  y1 ~~ y1
  y2 ~~ x2
  y1 ~~ x1
'
eu_clpm <- sem(CLPM_u, data = (gss_short1),missing="fiml",fixed.x=F)

auth_results <- rbind(parameterestimates(lr_clpm,level = .90),parameterestimates(eu_clpm,level = .90)) %>% filter(lhs %in% c("x2","x3","y2","y3") & rhs %in% c("x2","x1","y2","y1") & op == "~") %>% slice(c(2,3,6,7,10,11))
auth_results$lhs <- c("DV: Policy Preference","DV: Personality")
auth_results$rhs <- c("IV: Personality","IV: Policy Preference")
auth_results$label <- c(rep("Ideology",4),rep("EU",2))
auth_results$lhs <- factor(auth_results$lhs,levels=c("DV: Policy Preference","DV: Personality"))
auth_results$wave <- c(rep("Wave 2 DV~Wave 1 IV",2),rep("Wave 3 DV~Wave 2 IV",2),rep("Wave 2 DV~Wave 1 IV",2))

ggplot(auth_results,aes(x=label,y=est,shape=wave))+geom_pointrange(aes(ymin=ci.lower,ymax=ci.upper),position=position_dodge(width=.5))+facet_wrap(lhs~rhs)+geom_hline(yintercept = 0,lty=2)+theme_bw(base_size = 15)+labs(shape="Model")+ylab("Standardized Coefficient")+xlab("Policy Variable") + theme(legend.position="bottom")
ggsave("Results/xlaggedplotgesis_extraversion.pdf",width=8,height=6)

#SI 2.9 - Exploratory results for Agreeableness--------------------------
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1
  x3 + y3 ~ x2 + y2
  x3 ~~ x3
  x2 ~~ x2
  x1 ~~ x1
  y3 ~~ y3
  y2 ~~ y2
  y1 ~~ y1
  y3 ~~ x3
  y2 ~~ x2
  y1 ~~ x1
'

#Left-right
xvars <- colnames(gesis)[grep("lr", x = colnames(gesis))][c(1,2,3)]
yvars <- colnames(gesis)[grep("agre", x = colnames(gesis))][c(3,6,9)]
gesis %>% select(xvars, yvars) -> gss_short1
names(gss_short1)[1:3] <- paste0("x", 1:3)
names(gss_short1)[4:6] <- paste0("y", 1:3)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$x3<- scale(gss_short1$x3)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$y3<- scale(gss_short1$y3)

lr_clpm <- sem(CLPM_u, data = (gss_short1),missing="ML",fixed.x=F)

#EU
xvars <- colnames(gesis)[grep("antiEU", x = colnames(gesis))][c(1,2)]
yvars <- colnames(gesis)[grep("agre", x = colnames(gesis))][c(6,9)]
gesis %>% select(xvars, yvars) -> gss_short1
names(gss_short1)[1:2] <- paste0("x", 1:2)
names(gss_short1)[3:4] <- paste0("y", 1:2)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)

CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1
  x2 ~~ x2
  x1 ~~ x1
  y2 ~~ y2
  y1 ~~ y1
  y2 ~~ x2
  y1 ~~ x1
'
eu_clpm <- sem(CLPM_u, data = (gss_short1),missing="fiml",fixed.x=F)

auth_results <- rbind(parameterestimates(lr_clpm,level = .90),parameterestimates(eu_clpm,level = .90)) %>% filter(lhs %in% c("x2","x3","y2","y3") & rhs %in% c("x2","x1","y2","y1") & op == "~") %>% slice(c(2,3,6,7,10,11))
auth_results$lhs <- c("DV: Policy Preference","DV: Personality")
auth_results$rhs <- c("IV: Personality","IV: Policy Preference")
auth_results$label <- c(rep("Ideology",4),rep("EU",2))
auth_results$lhs <- factor(auth_results$lhs,levels=c("DV: Policy Preference","DV: Personality"))
auth_results$wave <- c(rep("Wave 2 DV~Wave 1 IV",2),rep("Wave 3 DV~Wave 2 IV",2),rep("Wave 2 DV~Wave 1 IV",2))

ggplot(auth_results,aes(x=label,y=est,shape=wave))+geom_pointrange(aes(ymin=ci.lower,ymax=ci.upper),position=position_dodge(width=.5))+facet_wrap(lhs~rhs)+geom_hline(yintercept = 0,lty=2)+theme_bw(base_size = 15)+labs(shape="Model")+ylab("Standardized Coefficient")+xlab("Policy Variable") + theme(legend.position="bottom")
ggsave("Results/xlaggedplotgesis_agreeableness.pdf",width=8,height=6)

#SI 2.10 -  Main results including demographic control----------------
#set model
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1 + variable + income4_1 + income4_2 + income4_3 + income4_4 + education_2 + education_3 + education_4 + education_5 + education_6 + age + female
  x3 + y3 ~ x2 + y2 + variable + income4_1 + income4_2 + income4_3 + income4_4 + education_2 + education_3 + education_4 + education_5 + education_6 + age + female
  x3 ~~ x3
  x2 ~~ x2
  x1 ~~ x1
  y3 ~~ y3
  y2 ~~ y2
  y1 ~~ y1
  y3 ~~ x3
  y2 ~~ x2
  y1 ~~ x1
  y1 ~ variable
'
#select data and prepare variables for left-right self-placement model
xvars <- colnames(all)[grep("lr", x = colnames(all))][c(1,2,3)]
yvars <- colnames(all)[grep("p", x = colnames(all))][c(8,9,10)]
all %>% select(xvars, yvars,z000001a,female,starts_with("education"),starts_with("income"),age,female,variable) -> gss_short1
names(gss_short1)[1:3] <- paste0("x", 1:3)
names(gss_short1)[4:6] <- paste0("y", 1:3)
gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$x3<- scale(gss_short1$x3)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$y3<- scale(gss_short1$y3)
#run model
lr_clpm <- sem(CLPM_u, data = (gss_short1),missing="ML",fixed.x=F, cluster="z000001a")

#set model for EU
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1 + variable + income4_1 + income4_2 + income4_3 + income4_4 + + education_2 + education_3 + education_4 + education_5 + education_6 + age + female
  x2 ~~ x2
  x1 ~~ x1
  y2 ~~ y2
  y1 ~~ y1
  y2 ~~ x2
  y1 ~~ x1
  y1 ~ variable
'

#select data and prepare variables for EU model
xvars <- colnames(all)[grep("antiEU", x = colnames(all))][c(1,2)]
yvars <- colnames(all)[grep("p", x = colnames(all))][c(9,10)]
all %>% select(xvars, yvars,z000001a,female,starts_with("education"),income,age,starts_with("income"),variable,age,female) -> gss_short1
names(gss_short1)[1:2] <- paste0("x", 1:2)
names(gss_short1)[3:4] <- paste0("y", 1:2)
gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
eu_clpm <- sem(CLPM_u, data = (gss_short1),missing="fiml",fixed.x=F,cluster = "z000001a")

#create figure
auth_results <- rbind(parameterestimates(lr_clpm,level = .90),parameterestimates(eu_clpm,level = .90)) %>% filter(lhs %in% c("x2","x3","y2","y3") & rhs %in% c("x2","x1","y2","y1") & op == "~") %>% slice(c(2,3,6,7,10,11))
auth_results$lhs <- c("DV: Policy Preference","DV: Personality")
auth_results$rhs <- c("IV: Personality","IV: Policy Preference")
auth_results$label <- c(rep("Ideology",4),rep("EU",2))
auth_results$lhs <- factor(auth_results$lhs,levels=c("DV: Policy Preference","DV: Personality"))
auth_results$wave <- c(rep("Wave 2 DV~Wave 1 IV",2),rep("Wave 3 DV~Wave 2 IV",2),rep("Wave 2 DV~Wave 1 IV",2))
ggplot(auth_results,aes(x=label,y=est,shape=wave))+geom_pointrange(aes(ymin=ci.lower,ymax=ci.upper),position=position_dodge(width=.5))+facet_wrap(lhs~rhs)+geom_hline(yintercept = 0,lty=2)+theme_bw(base_size = 15)+labs(shape="Model")+ylab("Standardized Coefficient")+xlab("Policy Variable") + theme(legend.position="bottom")
ggsave("Results/xlaggedplotgesiswithcontrols.pdf",width=8,height=6)

#SI 9.1 - Results conditional upon age---------

#left-righti
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1 + variable
  x3 + y3 ~ x2 + y2 + variable
  x3 ~~ x3
  x2 ~~ x2
  x1 ~~ x1
  y3 ~~ y3
  y2 ~~ y2
  y1 ~~ y1
  y3 ~~ x3
  y2 ~~ x2
  y1 ~~ x1
  y1 ~ variable
'

xvars <- colnames(all)[grep("lr", x = colnames(all))][c(1,2,3)]
yvars <- colnames(all)[grep("p", x = colnames(all))][c(8,9,10)]
all %>% select(xvars, yvars,id,female,starts_with("education"),starts_with("income"),age,variable) -> gss_short1
names(gss_short1)[1:3] <- paste0("x", 1:3)
names(gss_short1)[4:6] <- paste0("y", 1:3)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$x3<- scale(gss_short1$x3)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$y3<- scale(gss_short1$y3)

gss_short1$age1_median<-ifelse(gss_short1$age>median(gss_short1$age, na.rm=T), 1,0)
lr_clpm_young <- sem(CLPM_u, data = subset(gss_short1, age1_median==0),missing="fiml",fixed.x=F,cluster = "id")
lr_clpm_old <- sem(CLPM_u, data = subset(gss_short1, age1_median==1),missing="fiml",fixed.x=F,cluster = "id")
unique(gss_short1$id)

auth_results <-
  gtools::smartbind(
    parameterestimates(lr_clpm_young,level = .90)[c(2,14,28,40),],
    parameterestimates(lr_clpm_old,level = .90)[c(2,14,28,40),])

auth_results <- rbind(parameterestimates(lr_clpm_young,level = .90),parameterestimates(lr_clpm_old,level = .90)) %>% filter(lhs %in% c("x2","x3","y2","y3") & rhs %in% c("x2","x1","y2","y1") & op == "~") %>% slice(c(2,3,6,7,10,11,14,15))

auth_results$lhs <- c("DV: Policy Preference","DV: Personality")
auth_results$rhs <- c("IV: Personality","IV: Policy Preference")
auth_results$label <- c(rep("Ideology",4))
auth_results$lhs <- factor(auth_results$lhs,levels=c("DV: Policy Preference","DV: Personality"))
auth_results$wave <- c(rep("Wave 2 DV~Wave 1 IV",2),rep("Wave 3 DV~Wave 2 IV",2))
auth_results$age <- c(rep("Young",4),rep("Old",4))
auth_results$age <- relevel(factor(auth_results$age),ref="Young")

plot_lr <-ggplot(auth_results,aes(x=label,y=est,shape=wave, color=age))+geom_pointrange(aes(ymin=ci.lower,ymax=ci.upper),position=position_dodge(width=.5))+facet_wrap(lhs~rhs)+geom_hline(yintercept = 0,lty=2)+theme_bw()+labs(shape="Model")+ylab("Standardized Coefficient")+theme(axis.text.x = element_text(angle = -45,hjust=0))+xlab("Policy Variable") + theme(legend.position="none") + ggtitle("A. Ideology")


#EU
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1 + variable
  x2 ~~ x2
  x1 ~~ x1
  y2 ~~ y2
  y1 ~~ y1
  y2 ~~ x2
  y1 ~~ x1
  y1 ~ variable
'
xvars <- colnames(all)[grep("antiEU", x = colnames(all))][c(1,2)]
yvars <- colnames(all)[grep("p", x = colnames(all))][c(9,10)]
all %>% select(xvars, yvars,id,female,starts_with("education"),income,age,variable,starts_with("income")) -> gss_short1
names(gss_short1)[1:2] <- paste0("x", 1:2)
names(gss_short1)[3:4] <- paste0("y", 1:2)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$age1_median<-ifelse(gss_short1$age>median(gss_short1$age, na.rm=T), 1,0)

eu_clpm_young <- sem(CLPM_u, data = subset(gss_short1, age1_median==0),missing="fiml",fixed.x=F)
eu_clpm_old <- sem(CLPM_u, data = subset(gss_short1, age1_median==1),missing="fiml",fixed.x=F)

auth_results <-
  gtools::smartbind(
    parameterestimates(eu_clpm_young,level = .90)[c(2,14),],
    parameterestimates(eu_clpm_old,level = .90)[c(2,14),])

auth_results <- rbind(parameterestimates(eu_clpm_young,level = .90),parameterestimates(eu_clpm_old,level = .90)) %>% filter(lhs %in% c("x2","x3","y2","y3") & rhs %in% c("x2","x1","y2","y1") & op == "~") %>% slice(c(2,3,6,7,10,11,14,15))

auth_results$lhs <- c("DV: Policy Preference","DV: Personality")
auth_results$rhs <- c("IV: Personality","IV: Policy Preference")
auth_results$label <- c(rep("EU",4))
auth_results$lhs <- factor(auth_results$lhs,levels=c("DV: Policy Preference","DV: Personality"))
auth_results$wave <- c(rep("Wave 2 DV~Wave 1 IV",2),rep("Wave 3 DV~Wave 2 IV",2))
auth_results$age <- c(rep("Young",2),rep("Old",2))
auth_results$age <- relevel(factor(auth_results$age),ref="Young")

plot_eu <- ggplot(auth_results,aes(x=label,y=est,shape=wave, color=age))+geom_pointrange(aes(ymin=ci.lower,ymax=ci.upper),position=position_dodge(width=.5))+facet_wrap(lhs~rhs)+geom_hline(yintercept = 0,lty=2)+theme_bw()+labs(shape="Model")+ylab("Standardized Coefficient")+theme(axis.text.x = element_text(angle = -45,hjust=0))+xlab("Policy Variable") + theme(legend.position="none") + ggtitle("B. EU")

#create plot
comb_plot<-grid.arrange(plot_lr, plot_eu, nrow = 2)
ggsave(comb_plot, file="Results/GESIS_conditional_age.pdf")

#SI 9.2 - Results conditional upon education---------

#left-right
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1 + variable
  x3 + y3 ~ x2 + y2 + variable
  x3 ~~ x3
  x2 ~~ x2
  x1 ~~ x1
  y3 ~~ y3
  y2 ~~ y2
  y1 ~~ y1
  y3 ~~ x3
  y2 ~~ x2
  y1 ~~ x1
  y1 ~ variable
'
xvars <- colnames(all)[grep("lr", x = colnames(all))][c(1,2,3)]
yvars <- colnames(all)[grep("p", x = colnames(all))][c(8,9,10)]
all %>% select(xvars, yvars,id,female,starts_with("education"),starts_with("income"),age,variable) -> gss_short1
names(gss_short1)[1:3] <- paste0("x", 1:3)
names(gss_short1)[4:6] <- paste0("y", 1:3)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$x3<- scale(gss_short1$x3)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$y3<- scale(gss_short1$y3)

gss_short1$edu01 <-ifelse(gss_short1$education>3, 1,0)
lr_clpm_young <- sem(CLPM_u, data = subset(gss_short1, edu01==0),missing="fiml",fixed.x=F,cluster = "id")
lr_clpm_old <- sem(CLPM_u, data = subset(gss_short1, edu01==1),missing="fiml",fixed.x=F,cluster = "id")

auth_results <- rbind(parameterestimates(lr_clpm_young,level = .90),parameterestimates(lr_clpm_old,level = .90)) %>% filter(lhs %in% c("x2","x3","y2","y3") & rhs %in% c("x2","x1","y2","y1") & op == "~") %>% slice(c(2,3,6,7,10,11,14,15))

auth_results$lhs <- c("DV: Policy Preference","DV: Personality")
auth_results$rhs <- c("IV: Personality","IV: Policy Preference")
auth_results$label <- c(rep("Ideology",4))
auth_results$lhs <- factor(auth_results$lhs,levels=c("DV: Policy Preference","DV: Personality"))
auth_results$wave <- c(rep("Wave 2 DV~Wave 1 IV",2),rep("Wave 3 DV~Wave 2 IV",2))
auth_results$age <- c(rep("Low",4),rep("High",4))
auth_results$age <- relevel(factor(auth_results$age),ref="Low")

plot_lr <-ggplot(auth_results,aes(x=label,y=est,shape=wave, color=age))+geom_pointrange(aes(ymin=ci.lower,ymax=ci.upper),position=position_dodge(width=.5))+facet_wrap(lhs~rhs)+geom_hline(yintercept = 0,lty=2)+theme_bw()+labs(shape="Model")+ylab("Standardized Coefficient")+theme(axis.text.x = element_text(angle = -45,hjust=0))+xlab("Policy Variable") + theme(legend.position="none") + ggtitle("A. Ideology")

#EU
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1 + variable
  x2 ~~ x2
  x1 ~~ x1
  y2 ~~ y2
  y1 ~~ y1
  y2 ~~ x2
  y1 ~~ x1
  y1 ~ variable
'
xvars <- colnames(all)[grep("antiEU", x = colnames(all))][c(1,2)]
yvars <- colnames(all)[grep("p", x = colnames(all))][c(9,10)]
all %>% select(xvars, yvars,id,female,starts_with("education"),income,age,variable,starts_with("income")) -> gss_short1
names(gss_short1)[1:2] <- paste0("x", 1:2)
names(gss_short1)[3:4] <- paste0("y", 1:2)

gss_short1$x1<- scale(gss_short1$x1)
gss_short1$x2<- scale(gss_short1$x2)
gss_short1$y1<- scale(gss_short1$y1)
gss_short1$y2<- scale(gss_short1$y2)
gss_short1$edu01 <-ifelse(gss_short1$education>3, 1,0)

eu_clpm_young <- sem(CLPM_u, data = subset(gss_short1, edu01==0),missing="fiml",fixed.x=F)
eu_clpm_old <- sem(CLPM_u, data = subset(gss_short1, edu01==1),missing="fiml",fixed.x=F)

auth_results <- rbind(parameterestimates(eu_clpm_young,level = .90),parameterestimates(eu_clpm_old,level = .90)) %>% filter(lhs %in% c("x2","x3","y2","y3") & rhs %in% c("x2","x1","y2","y1") & op == "~") %>% slice(c(2,3,6,7,10,11,14,15))

auth_results$lhs <- c("DV: Policy Preference","DV: Personality")
auth_results$rhs <- c("IV: Personality","IV: Policy Preference")
auth_results$label <- c(rep("EU",4))
auth_results$lhs <- factor(auth_results$lhs,levels=c("DV: Policy Preference","DV: Personality"))
auth_results$wave <- c(rep("Wave 2 DV~Wave 1 IV",2),rep("Wave 3 DV~Wave 2 IV",2))
auth_results$age <- c(rep("Low",2),rep("High",2))
auth_results$age <- relevel(factor(auth_results$age),ref="Low")

plot_eu <- ggplot(auth_results,aes(x=label,y=est,shape=wave, color=age))+geom_pointrange(aes(ymin=ci.lower,ymax=ci.upper),position=position_dodge(width=.5))+facet_wrap(lhs~rhs)+geom_hline(yintercept = 0,lty=2)+theme_bw()+labs(shape="Model")+ylab("Standardized Coefficient")+theme(axis.text.x = element_text(angle = -45,hjust=0))+xlab("Policy Variable") + theme(legend.position="none") + ggtitle("B. EU")

#create plot
comb_plot<-grid.arrange(plot_lr, plot_eu, nrow = 2)
ggsave(comb_plot, file="Results/GESIS_conditional_edu.pdf")


